source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

# import data ---------
library(GEOquery)
getGEOSuppFiles('GSE212092', makeDirectory = F, fetch_files = F)

platf_meta <- getGEO('GSE212092')

platf_meta$`GSE212092-GPL18573_series_matrix.txt.gz`$`cell type:ch1` |> unique()
platf_meta$`GSE212092-GPL24676_series_matrix.txt.gz`$`tissue:ch1` |> unique()

pbmc_path <- list.files('mission/SLE_TRPM2_MfMo/GSE212092_LPS/', 'PBMC',
                        full.names = T)

bm_path <- list.files('mission/SLE_TRPM2_MfMo/GSE212092_LPS/', '_BM',
                       full.names = T)

pbmc_mex <- pbmc_path |>
  read_geo_supp('Subject\\d_[:alnum:]+')

sobj <- bm_path |>
  read_geo_supp('Subject\\d_[:alnum:]+')

pbmc_mex <- pbmc_mex |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

pbmc_mex |> VlnPlot('mito_ratio', pt.size = 0)

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- pbmc_mex |> filter(mito_ratio < 10)

sobj <- sobj |> filter(mito_ratio < 10)

sobj <- sobj |>
  quick_process_seurat()

bm_markers <- c('HBD','SPINK2','PF4','CD68','MAFB')

sobj |>
  DotPlot(c(pbmc_markers, bm_markers), cols = 'RdYlBu', cluster.idents = T) +
  RotatedAxis()

pbmc_cluster_excel <-
  read_tsv('mission/SLE_TRPM2_MfMo/lps_pbmc_cluster_excel.tsv')

bm_cluster_excel <-
  read_tsv('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_bm_cluster_excel.tsv')

sobj <- sobj |>
  mutate(cluster_id = as.double(seurat_clusters)) |>
  left_join(pbmc_cluster_excel)

sobj <- sobj |>
  mutate(seurat_clusters = as.double(seurat_clusters)) |>
  left_join(bm_cluster_excel)

sobj <- sobj |>
  filter(manual_main != 'uk') |>
  mutate(manual_main = case_match(manual_main,
                                  'B' ~ 'B cell',
                                  'mast' ~ 'Mast cell',
                                  'NK' ~ 'NK cell',
                                  'PC' ~ 'Plasma cell',
                                  'myeloid' ~ 'Myeloid cell',
                                  .default = manual_main))

sobj <- sobj |> filter(manual_main != 'doublet')

sobj |>
  DimPlot(group.by = 'manual_main', cols = 'Paired')

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'Day\\d|4h') |>
           case_match('Day0' ~ '0 hpi',
                      '4h' ~ '4 hpi',
                      'Day7' ~ '7 dpi'))

sobj <- sobj |>
  mutate(manual_fine = case_when(cluster_id %in% c(15,20,18) ~ 'DC',
                                 manual_main == 'Myeloid cell' ~ 'Monocyte',
                                 .default = manual_main))

# big umap -----------
sobj |>
  DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('PBMC') +
  theme_jpub(theme_classic)

publish_pdf('LPS.PBMC.umap.pdf', width = 70)

sobj |>
  filter(manual_main != 'pDC') |>
  DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('Bone marrow') +
  theme_jpub(theme_classic)

publish_pdf('LPS.BM.umap.pdf', width = 70)

# TRPM2 expr overview ------------
sobj |>
  DotPlot2d('TRPM2', group, manual_fine) +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Time')

trpm2_pbmc_time <- last_plot()

trpm2_pbmc_time$data |>
  mutate(group.y = fct_reorder(group.y, avg.exp, max)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Cell type', y = 'Time post LPS injection',
       title = 'TRPM2 expression in PBMC') +
  RotatedAxis() +
  theme_jpub()

publish_source_plot('LPS.pbmc.TRPM2.dotplot', width = 60)

trpm2_bm_time <- last_plot()

trpm2_bm_time$data |>
  mutate(group.y = fct_reorder(group.y, avg.exp, max)) |>
  filter(group.y != 'pDC') |>
  BubblePlot(d2 = T) +
  labs(x = 'Cell type', y = 'Time post LPS injection',
       title = 'TRPM2 expression in BM') +
  RotatedAxis() +
  theme_jpub()

publish_source_plot('LPS.BM.TRPM2.dotplot', width = 60)

sobj_mo$trpm2_type |> table()

sobj_mo |>
  DotPlot(c('TRPM2','CD14'), cols = 'RdBu')

sobj$manual_main <- sobj$manual_fine

m2h_bc <- sobj_mo |>
  filter(manual_fine != 'ncMo') |>
  colnames()

sobj <- sobj |>
  mutate(manual_fine = case_when(.cell %in% m2h_bc ~ 'CD14+ Monocyte',
                                 manual_main == 'Monocyte' ~ 'CD16+ Monocyte',
                                 .default = manual_main))

# full rds ----------------
sobj |>
  write_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_pbmc.rds')

sobj |>
  write_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_BM.rds')

sobj <-
  read_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_pbmc.rds')

sobj <-
  read_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_BM.rds')

sobj <- sobj |>
  filter(manual_main != 'pDC')

# monocyte subset -----------
sobj_mo <- sobj |>
  filter(manual_fine == 'Monocyte')

sobj_mo <- sobj |>
  filter(manual_main == 'Mo/DC')

sobj_mo <- sobj_mo |>
  quick_process_seurat(skip_norm = T)

modc_marker <- list('cMono'=c('CD14','S100A8','S100A9'),
                    'ncMono' = c('FCGR3A','CDKN1C'),
                    'Macrophage'=c('CD68','MAFB','MARCO','CD163','MRC1'),
                    'cDC'=c('FCER1A','CLEC9A','XCR1','CD1C','CLEC10A'))

sobj_mo |>
  DotPlot(c('TRPM2', modc_marker), cols = 'RdYlBu',
          cluster.idents = T, dot.scale = 3) +
  labs(x = 'Gene', y = 'Cell cluster') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  RotatedAxis()

g1 <- last_plot()

g1$data |> BubblePlot() +
  labs(x = 'Gene', y = 'Cell cluster') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('lps.BM.momfdc.marker.dotplot', width = 80)

sobj_mo <- sobj_mo |>
  filter(seurat_clusters != 14)

sobj_mo <- sobj_mo |>
  mutate(manual_fine = case_when(seurat_clusters %in% c(10,13) ~ 'ncMo',
                                 seurat_clusters == 7 ~ 'cDC',
                                 seurat_clusters %in% c(8,1) ~ 'Mo-MF',
                                 .default = 'cMo'))

sobj_mo <- sobj_mo |>
  filter(str_detect(manual_fine, 'Mo'), nCount_RNA < 40000)

sobj_mo <- sobj_mo |>
  quick_process_seurat(skip_norm = T)

## BM subtype -------------
sobj_mo <- sobj_mo |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(1,2,9),
                             'TRPM2-hi Mono', 'TRPM2-lo Mono'))

## mono rds -----------
sobj_mo <- sobj_mo |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(5),
                             'TRPM2-lo Mono', 'TRPM2-hi Mono'))

sobj_mo |>
  write_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_pbmc_mono.rds')

sobj_mo |>
  write_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_BM_mono.rds')

sobj_mo <-
  read_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_pbmc_mono.rds')

sobj_mo <-
  read_rds('mission/SLE_TRPM2_MfMo/GSE212092_LPS/lps_BM_mono.rds')

sobj <- sobj |>
  mutate(manual_main = case_when(.cell %in% colnames(sobj_mo) ~ 'Monocyte',
                                 manual_main == 'Mo/DC' ~ 'DC',
                                 .default = manual_main))

sobj_mo |>
  DimPlot(cols = 'Paired', label = T, label.box = T, repel = T,
          label.size = 2) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_pdf('LPS.PBMC.mono.umap.pdf', width = 60)

publish_pdf('LPS.BM.mono.umap.pdf', width = 60)

sobj_mo |>
  DimPlot(group.by = 'manual_fine') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'Monocyte subsets') +
  theme_jpub

publish_pdf('LPS.PBMC.mono.m2subset.umap.pdf', width = 60)

sobj_mo |>
  DimPlot(group.by = 'trpm2_type') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'Monocyte/macrophage subsets') +
  theme_jpub

publish_pdf('LPS.BM.mono.m2subset.umap.pdf', width = 60)

sobj_mo |>
  FeaturePlot('TRPM2', order = T, split.by = 'group', pt.size = .05,
              cols = c('lightgrey','red')) &
  theme_classic(base_size = 6, base_family = 'ArialMT') &
  NoLegend()

publish_pdf('LPS.PBMC.mono.TRPM2.featureplot.pdf', width = 100)

publish_pdf('LPS.BM.mono.TRPM2.featureplot.pdf', width = 120)

sobj_mo |>
  DotPlot('TRPM2', cols = 'RdYlBu', dot.scale = 2) +
  labs(x = 'Gene', y = 'Monocyte clusters') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('LPS.PBMC.mono.TRPM2.dotplot')

publish_source_plot('LPS.BM.mono.TRPM2.dotplot')

### TNF source? ------------
sobj |>
  DotPlot2d('TNF', group, manual_fine) +
  labs(x = 'Time', y = 'Cell type', title = 'TNF in PBMC') +
  theme_jpub()

publish_source_plot('lps.pbmc.timepoint.tnf.dotplot')

sobj |>
  filter(manual_main != 'pDC') |>
  DotPlot2d('TNF', group, manual_fine) +
  labs(x = 'Time', y = 'Cell type', title = 'TNF in BM') +
  theme_jpub()

publish_source_plot('lps.bm.timepoint.tnf.dotplot')

tnf_logfc <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi', features = 'TNF')

tnf_logfc |> DT::datatable()

sobj_mo |>
  DotPlot('TNF', group.by = 'orig.ident')

mo_tnf <- last_plot() |>
  pluck('data') 

mo_tnf |>
  mutate(time = case_when(str_detect(id, '4h') ~ '4 hpi',
                          str_detect(id, 'y0') ~ '0 hpi',
                          .default = '7 dpi')) |>
  ggplot(aes(time, avg.exp)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(method = 't.test', ref.group = '0 hpi')

sobj_mo |>
  bill.violin('TNF', group.by = group)

### CCR2 profile -------------
sobj |>
  DotPlot2d('CCR2', group, manual_main) +
  labs(x = 'Time', y = 'Cell type', title = 'CCR2 in BM')

sobj |>
  filter(manual_main != 'pDC') |>
  DotPlot2d('CCR2', group, manual_main) +
  labs(x = 'Time', y = 'Cell type', title = 'CCR2 in BM')

### TRP family -----------
trp_family <- rownames(sobj) |> str_subset('^TRP(C|M|V)\\d$') |>
  str_sort()

sobj |>
  filter(group == '0 hpi') |>
  DotPlot(trp_family, group.by = 'manual_fine', cols = 'RdBu', dot.scale = 3) +
  labs(x = 'Gene', y = 'Cell type', title = 'BM baseline') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('lps.pbmc.0hpi.trp.family.dotplot', width = 80)

publish_source_plot('lps.bm.0hpi.trp.family.dotplot', width = 80)

sobj |>
  filter(group == '4 hpi') |>
  DotPlot(trp_family, group.by = 'manual_fine', cols = 'RdBu', dot.scale = 3) +
  labs(x = 'Gene', y = 'Cell type', title = 'BM 4 hpi LPS') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('lps.pbmc.4hpi.trp.family.dotplot', width = 80)

publish_source_plot('lps.bm.4hpi.trp.family.dotplot', width = 80)

sobj |>
  filter(group == '7 dpi') |>
  DotPlot(trp_family, group.by = 'manual_fine', cols = 'RdBu', dot.scale = 3) +
  labs(x = 'Gene', y = 'Cell type', title = 'BM 7 dpi LPS') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('lps.pbmc.7dpi.trp.family.dotplot', width = 80)

publish_source_plot('lps.bm.7dpi.trp.family.dotplot', width = 80)

pbmc_subtype_lps_4hpi_deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi')

pbmc_subtype_lps_4hpi_trp <- pbmc_subtype_lps_4hpi_deg |>
  mutate(cluster = as.character(cluster)) |>
  filter(gene %in% trp_family) 

pbmc_subtype_lps_4hpi_trp |>
  ggplot(aes(gene, cluster, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape = 21) +
  scale_size() +
  theme_jpub() +
  labs(y = 'Cell type', title = 'log2FC in PBMC 4 hpi LPS') +
  RotatedAxis() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(pbmc_subtype_lps_4hpi_trp$avg_log2FC))

publish_source_plot('lps.pbmc.4hpi.trp.family.logfc', width = 70)

pbmc_subtype_lps_7dpi_deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'group',
                       ident.1 = '7 dpi', ident.2 = '0 hpi')

pbmc_subtype_lps_7dpi_trp <- pbmc_subtype_lps_7dpi_deg |>
  mutate(cluster = as.character(cluster)) |>
  filter(gene %in% trp_family) 

pbmc_subtype_lps_7dpi_trp |>
  ggplot(aes(gene, cluster, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape = 21) +
  scale_size() +
  theme_jpub() +
  labs(y = 'Cell type', title = 'log2FC in PBMC 7 dpi LPS') +
  RotatedAxis() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(pbmc_subtype_lps_7dpi_trp$avg_log2FC))

publish_source_plot('lps.pbmc.7dpi.trp.family.logfc', width = 70)

## TRPM2 fold change between time points -----------
acute_mo_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'trpm2_type', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi')

acute_mo_deg |>
  filter(gene == 'TRPM2')

sobj_mo |>
  filter(trpm2_type == 'TRPM2-hi Mo-Mac') |>
  bill.violin('TRPM2', group.by = group) +
  theme_bw(base_size = 6, base_family = 'ArialMT')

g1 <- last_plot()

g3 <- g1 + NoLegend() +
  labs(subtitle = 'TRPM2-high monocyte', y = 'Normalized expression',
       title = 'TRPM2 expression in PBMC monocytes after LPS injection') +
  expand_limits(y = 2.9)

g2 <- sobj_mo |>
  filter(manual_fine == 'TRPM2-lo Mono') |>
  bill.violin('TRPM2', group.by = group) +
  theme_bw(base_size = 6, base_family = 'ArialMT')

g4 <- g2 + NoLegend() +
  labs(subtitle = 'TRPM2-low monocyte', y = 'Normalized expression') +
  expand_limits(y = 2.9)

g3 + g4

publish_pdf('LPS.PBMC.mo.TRPM2.time.violin.pdf', width = 70)

sobj_mo |>
  bill.violin('TRPM2', group.by = group) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  NoLegend() +
  labs(subtitle = 'Monocyte', y = 'Normalized expression')

publish_pdf('LPS.PBMC.totalmono.TRPM2.time.violin.pdf', width = 40)

publish_pdf('LPS.BM.totalmono.TRPM2.time.violin.pdf', width = 40)

## 4h vs 0h TRPM2 logfc in leiden clusters ---------
pbmc_leiden_m2_4v0 <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi', features = 'TRPM2')

pbmc_leiden_m2_4v0 |>
  ggplot(aes(cluster, fill = -log10(p_val_adj), y = avg_log2FC)) +
  geom_col() +
  theme_bw() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(x = 'Monocyte clusters', title = 'TRPM2 logfc in LPS PBMC')

bm_leiden_m2_4v0 <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi', features = 'TRPM2')

bm_leiden_m2_4v0 |>
  mutate(cluster = fct_reorder(cluster, avg_log2FC)) |>
  ggplot(aes(cluster, fill = -log10(p_val_adj), y = avg_log2FC)) +
  geom_col() +
  theme_bw() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(x = 'Monocyte clusters', title = 'TRPM2 logfc in LPS BM')

## M2-hi subset fraction ------------
count_indiv <- sobj |>
  summarise(sum = n(), .by = orig.ident)

count_indiv

frac_m2h_in_mono <- sobj_mo |>
  calc_frac_conf_on_grouped_count(orig.ident, trpm2_type) |>
  mutate(group = str_extract(orig.ident, '4h|Day.') |>
           case_match('Day0' ~ '0 hpi',
                      '4h' ~ '4 hpi',
                      'Day7' ~ '7 dpi')) |>
  filter(str_detect(trpm2_type, 'hi'))

frac_m2h_in_mono |>
  ggplot(aes(group, 100*fraction, color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('0 hpi', '4 hpi'), c('0 hpi', '7 dpi')),
                     method = 't.test', size = 1.5) +
  labs(x = 'Time after LPS injection', y = '% in total monocyte',
       title = 'Proportion of TRPM2-high monocyte') +
  NoLegend() +
  theme_bw(base_size = 6, base_family = 'ArialMT')

wide_frac <- frac_m2h_in_mono |>
  mutate(indiv = str_extract(orig.ident, 'Subject.')) |>
  pivot_wider(names_from = group, values_from = fraction, id_cols = indiv)

t.test(wide_frac$`0 hpi`, wide_frac$`4 hpi`, paired = T)

publish_source_plot('LPS_PBMC_M2hi_mo_frac_in_totalmono')

publish_source_plot('LPS_BM_M2hi_mo_frac_in_totalmono')

frac_m2h_in_mono <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS_PBMC_M2hi_mo_frac_in_totalmono.csv')

frac_m2h_in_mono |>
  ggplot(aes(group, 100*(1-fraction), color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('0 hpi', '4 hpi'), c('0 hpi', '7 dpi')),
                     method = 't.test', size = 1.5) +
  labs(x = 'Time after LPS injection', y = '% in total monocyte',
       title = 'Proportion of TRPM2-low monocyte') +
  NoLegend() +
  theme_bw(base_size = 6, base_family = 'ArialMT')

publish_source_plot('LPS_PBMC_M2lo_mo_frac_in_totalmono')

publish_source_plot('LPS_BM_M2lo_mo_frac_in_totalmono')

frac_m2h_in_mono |>
  left_join(count_indiv) |>
  mutate(frac_in_pmbc = n/sum) |>
  ggplot(aes(group, frac_in_pmbc)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('0 hpi', '4 hpi'), c('0 hpi', '7 dpi')),
                     method = 't.test') +
  labs(x = 'Time after LPS injection', y = '% in PBMC')

## inflammatory module ---------
infl_mod <- read_delim('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sobj_mo <- sobj_mo |>
  AddModuleScore(features = list(infl_mod$gene), name = 'inflam')

tag_gene_inflam <- sobj_mo |>
  DotPlot(c('TRPM2','HMGB2','CCR2','inflam1','CD14'), cols = 'RdYlBu') |>
  pluck('data') |>
  as_tibble() |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(Subtype = ifelse(TRPM2 < .1, 'TRPM2-lo Mono', 'TRPM2-hi Mono')) 

tag_gene_inflam |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = Subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.trpm2.inflam.correlation', width = 70)

publish_source_plot('LPS.BM.mo.trpm2.inflam.correlation', width = 70)

tag_gene_inflam |>
  ggplot(aes(CCR2, inflam1)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = Subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Average expression of CCR2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.CCR2.inflam.correlation', width = 70)

tag_gene_inflam |>
  ggplot(aes(HMGB2, inflam1)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Average expression of HMGB2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.HMGB2.inflam.correlation', width = 70)

tag_gene_inflam <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS.pbmc.mo.HMGB2.inflam.correlation.csv')

tag_gene_inflam |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Average expression of CCR2', x = 'Average expression of HMGB2') +
  theme_jpub()

publish_source_plot('LPS.pbmc.mo.TRPM2.CCR2.correlation', width = 70)

tag_gene_inflam <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS.BM.mo.CD14.TRPM2.correlation.csv')

tag_gene_inflam |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Average expression of CCR2', x = 'Average expression of HMGB2') +
  theme_jpub()

publish_source_plot('LPS.BM.mo.TRPM2.CCR2.correlation', width = 70)

### added scaled gene --------
scaled_gene_inflam <- tag_gene_inflam |>
  mutate(scaled_TRPM2 = scale(TRPM2)[,1],
         scaled_CCR2 = scale(CCR2)[,1],
         scaled_HMGB2 = scale(HMGB2)[,1])

scaled_gene_inflam |>
  ggplot(aes(scaled_TRPM2 + scaled_CCR2, inflam1)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Scaled expression of TRPM2+CCR2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.TRPM2_CCR2.inflam.correlation', width = 70)

scaled_gene_inflam |>
  ggplot(aes(scaled_TRPM2 + scaled_HMGB2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = Subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Scaled expression of TRPM2+HMGB2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.TRPM2_HMGB2.inflam.correlation', width = 70)

publish_source_plot('LPS.BM.mo.TRPM2_HMGB2.inflam.correlation', width = 70)

scaled_gene_inflam <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS.pbmc.mo.TRPM2_HMGB2.inflam.correlation.csv')

scaled_gene_inflam |>
  ggplot(aes(scaled_HMGB2 + scaled_CCR2, inflam1)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Scaled expression of HMGB2+CCR2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.HMGB2_CCR2.inflam.correlation', width = 70)

scaled_gene_inflam |>
  ggplot(aes(scaled_HMGB2 + scaled_CCR2, inflam1)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Scaled expression of HMGB2+CCR2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.HMGB2_CCR2.inflam.correlation', width = 70)

scaled_gene_inflam |>
  ggplot(aes(scaled_HMGB2 + scaled_CCR2 + scaled_TRPM2, inflam1)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'PBMC monocyte clusters after LPS administration',
       y = 'Inflammatory module score', x = 'Scaled expression of TRPM2+HMGB2+CCR2') +
  theme_jpub

publish_source_plot('LPS.pbmc.mo.TRPM2_HMGB2_CCR2.inflam.correlation', width = 70)

### CD14 x TRPM2 -----------
tag_gene_inflam |>
  ggplot(aes(TRPM2, CD14)) +
  geom_point(aes(color = Subtype)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'BM monocyte clusters after LPS administration') +
  theme_jpub()

publish_source_plot('LPS.pbmc.mo.CD14.TRPM2.correlation', width = 70)

publish_source_plot('LPS.BM.mo.CD14.TRPM2.correlation', width = 70)


## inflam violin ----------
pbmc_4v0_m2h_inf12 <- lps_4hpi_mo_deg |>
  filter(gene %in% pbmc_4v0_m2h_infgene, cluster == 'TRPM2-hi Mono') |>
  slice_min(p_val_adj, n = 12)

sobj_mo |>
  bill.violin(c('Hmgb2','Cd14','Ccr2', pbmc_4v0_m2h_inf12$gene),
              group.by = trpm2_type, facet.ncol = 5, pubsize = T) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in monocyte 4 hpi LPS')

publish_pdf('LPS.4hpi.m2hl.inflam.gene.violin.pdf', width = 90)

bm_4v0_m2h_inf12 <- lps_4hpi_mo_deg |>
  filter(gene %in% bm_4v0_m2h_infgene, cluster == 'cMo') |>
  slice_min(p_val_adj, n = 12)

sobj_mo |>
  bill.violin(c('Hmgb2','Cd14','Ccr2', bm_4v0_m2h_inf12$gene),
              group.by = trpm2_type, facet.ncol = 5, pubsize = T) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in monocyte 4 hpi LPS')

publish_pdf('LPS.BM.4hpi.m2hl.inflam.gene.violin.pdf', width = 90)

sobj_mo |>
  filter(str_detect(trpm2_type, 'hi')) |>
  bill.violin(c('Hmgb2','Cd14','Ccr2', pbmc_4v0_m2h_inf12$gene),
              group.by = group, facet.ncol = 5, pubsize = T) +
  labs(x = 'Time', y = 'Normalized expression', fill = 'Time',
       title = 'Inflammatory response genes in monocyte 4 hpi LPS')

publish_pdf('LPS.m2hi.time.inflam.gene.violin.pdf', width = 90)

sobj_mo |>
  filter(str_detect(trpm2_type, 'hi')) |>
  bill.violin(c('Hmgb2','Cd14','Ccr2', bm_4v0_m2h_inf12$gene),
              group.by = group, facet.ncol = 5, pubsize = T) +
  labs(x = 'Time', y = 'Normalized expression', fill = 'Time',
       title = 'Inflammatory response genes in monocyte 4 hpi LPS')

publish_pdf('LPS.BM.m2hi.time.inflam.gene.violin.pdf', width = 90)

## TAB-MAP3K7 program ------------
tab_complex <- c('MAP3K7', 'TAB1', 'TAB2', 'TAB3')

m2hvl_timesplit_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'group', group.by = 'manual_fine',
                       ident.1 = 'TRPM2-hi Mono')

m2hvl_timesplit_deg |>
  filter(gene %in% tab_complex, p_val_adj < .05)

m2hvl_deg <- sobj_mo |>
  FindMarkers(group.by = 'manual_fine', ident.1 = 'TRPM2-hi Mono')

m2hvl_deg |>
  as_tibble(rownames = 'gene') |>
  filter(gene %in% tab_complex, p_val_adj < .05)

sobj_mo |>
  bill.violin(tab_complex, group.by = manual_fine, facet.ncol = 2) +
  labs(x = 'Cell type', fill = 'Cell type', y = 'Normalized expression',
       title = 'TAK1 complex genes in monocyte') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_source_plot('LPS.pbmc.mono.TAB-MAP3K7.violin', width = 60)

sobj_mo |>
  filter(group == '4 hpi') |>
  mutate(sample = str_c(manual_fine, '_', orig.ident)) |>
  DotPlot(group.by = 'sample', tab_complex, cols = 'RdYlBu')

## GSEA 4hpi vs 0hpi ---------
library(clusterProfiler)

lps_4hpi_mo_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi')

m2h_4hpi_gsego <- acute_mo_deg |>
  filter(str_detect(cluster, 'hi'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'BP', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0, pvalueCutoff = 1)

m2h_4hpi_simgo <- m2h_4hpi_gsego |> simplify()

m2h_4hpi_gsego@result |>
  filter(NES > 0, qvalue < .05) |>
  slice_min(p.adjust, n = 20) |>
  slice_sample(n = 10) |>
  plot_enrichment(metric = NES) +
  theme_pubr(base_size = 6, base_family = 'ArialMT', legend = 'right') +
  labs(x = 'Normalized enrichment score', y = 'Pathway',
       title = 'GO BP pathway enrichment of TRPM2-hi mono:\n4 hpi LPS vs 0 hpi') +
  theme_jpub

publish_source_plot('LPS.4hpi.vs.0hpi.pbmc.m2h.gogse', width = 100)

pbmc_4v0_m2h_gogse <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS.4hpi.vs.0hpi.pbmc.m2h.gogse.csv')

pbmc_4v0_m2h_infgene <- pbmc_4v0_m2h_gogse |>
  filter(str_starts(Description, 'infla')) |>
  pull(core_enrichment) |>
  str_split_1('/')

publish_source_plot('LPS.4hpi.vs.0hpi.BM.m2h.gogse', width = 100)

bm_4v0_m2h_gogse <-
  read_csv('mission/SLE_TRPM2_MfMo/results/LPS.4hpi.vs.0hpi.BM.m2h.gogse.csv')

bm_4v0_m2h_infgene <- bm_4v0_m2h_gogse |>
  filter(str_starts(Description, 'infla')) |>
  pull(core_enrichment) |>
  str_split_1('/')

lps_4hpi_m2mo_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'trpm2_type', group.by = 'group',
                       ident.1 = '4 hpi', ident.2 = '0 hpi')

lps_4hpi_m2mo_deg |>
  filter(str_detect(gene, 'NFKB\\d'))

g1 <- sobj_mo |>
  filter(trpm2_type == 'TRPM2-hi Mono') |>
  bill.violin(c('NFKB1','NFKB2'), group.by = group) +
  labs(title = 'TRPM2-high monocyte') +
  RotatedAxis()

g2 <- sobj_mo |>
  filter(trpm2_type == 'TRPM2-lo Mono') |>
  bill.violin(c('NFKB1','NFKB2'), group.by = group) +
  labs(title = 'TRPM2-low monocyte') +
  RotatedAxis()

g1 + g2

### type 1 IFN ----------
library(enrichplot)

m2h_4hpi_gsego |>
  gseaplot2('GO:0034340',
            title = 'Response to type I interferon (NES=-0.90, FDR=0.68)')

m2h_4hpi_gsego@result |>
  filter(str_detect(Description, ' I interferon'))

t1ifn <- sobj_mo |> rownames() |> str_subset('^IFN(A|B|W|E|K)')

lps_4hpi_mo_deg |>
  filter(str_detect(gene, 'MX1|IRF7'))

## GSEA 4h m2h vs m2l -------------
m2hvl_bytime_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'group', group.by = 'manual_fine',
                       ident.1 = 'TRPM2-hi Mono')

m2hvl_4hpi_gsego <- m2hvl_bytime_deg |>
  filter(cluster == '4 hpi', p_val_adj < 0.05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL',
        eps = 0, pvalueCutoff = 1)

m2hvl_4hpi_gsego@result |>
  filter(ONTOLOGY == 'BP') |>
  as_tibble()
  plot_enrichment(metric = NES) +
  theme_pubr(base_size = 6, base_family = 'ArialMT', legend = 'right') +
  labs(x = 'Normalized enrichment score', y = 'Pathway',
       title = 'GO BP pathway enrichment of 4 hpi LPS:\nTRPM2-hi mono vs TRPM2-lo mono') +
  theme_jpub

publish_source_plot('LPS.4hpi.vs.0hpi.pbmc.m2h.gogse', width = 100)

### type 1 IFN ----------
library(enrichplot)

m2h_4hpi_gsego |>
  gseaplot2('GO:0034340',
            title = 'Response to type I interferon (NES=-0.90, FDR=0.68)')

m2h_4hpi_gsego@result |>
  filter(str_detect(Description, ' I interferon'))

resp2_t1ifn <- map_go_gene('GO:0034340')

sig_resp2_t1ifn <- m2hvl_bytime_deg |>
  filter(cluster == '4 hpi', gene %in% resp2_t1ifn$SYMBOL, p_val_adj < .05)

select <- dplyr::select

sobj_mo |>
  filter(group == '4 hpi') |>
  bill.violin(sig_resp2_t1ifn$gene, group.by = manual_fine, pubsize = T) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Response to type I IFN genes in monocyte 4 hpi LPS')

## GSEA m2h vs m2l ------------
m2hvl_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'group', group.by = 'manual_fine', ident.1 = 'TRPM2-hi Mono')

m2hvl_deg |>
  filter(gene == 'CCR2')

chemotax_gene <- search_go_term('^chemotaxis$') |>
  pull(GOID) |>
  map_go_gene()

m2hvl_deg |>
  filter(cluster == '0 hpi', gene %in% chemotax_gene$SYMBOL, avg_log2FC > 0) |>
  slice_min(p_val_adj, n = 20) |>
  mutate(
    p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
    gene = fct_reorder(gene, avg_log2FC)
  ) |>
  ggplot(aes(avg_log2FC, gene, fill = -log10(p_val_adj))) +
  geom_col() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  theme_jpub() +
  labs(title = 'Top enriched chemotaxis genes:\nTRPM2-high vs TRPM2-low monocytes')

publish_source_plot('LPS.PBMC.chemotaxis.gene.barplot')

## pseudotime analysis ---------
library(slingshot)

sobj_mo |> DimPlot()

sobj_mo |> DimPlot(group.by = 'manual_fine')

foo <- sobj_mo |>
  slingshot_seurat(lineage = 'Lineage1', cluster = 'manual_fine')

g1 <- last_plot()

g1 + theme_jpub(theme_classic)

publish_pdf('lps.bm.mo.pseudotime.umap.pdf', width = 70)

m2_percell <- foo |>
  get_abundance_sc_wide(c('TRPM2','CD14'))

foo |>
  select(Lineage1, manual_fine) |>
  bind_cols(m2_percell) |>
  ggplot(aes(Lineage1, TRPM2)) +
  geom_point(aes(y = TRPM2/15, color = manual_fine), size = AutoPointSize(m2_percell)) +
  geom_smooth() +
  labs(x = 'Pseudotime', color = 'Subtype',
       title = 'Expression of TRPM2 along infered trajectory') +
  theme_jpub()

publish_source_plot('lps.bm.mo.pseudotime.trpm2.infer')

Misc(foo)

sling_mo <- sobj_mo |>
  Embeddings(reduction = 'umap') |>
  slingshot(sobj_mo$manual_fine)

mo_pt <- sling_mo |>
  slingPseudotime()

mo_pt <- mo_pt |> as.data.frame()

mo_pt |> select('Lineage1') |> head()

sds <- as.SlingshotDataSet(osteo_ss)

sds |> glimpse()

sling_mo |> slingshot::slingCurves() |>
  map('s')

sling_mo@metadata$curves$Lineage1$s |> head()

sds@curves$Lineage1$s |>
  head()

sobj |>
  AddMetaData(osteo_pt) |>
  FeaturePlot('Lineage1', cols = c('lightgrey','red')) +
  ggtitle('Pseudotime: Lineage1 in Prrx1+ cells') +
  geom_path(data = sds@curves$Lineage1$s, aes(x = umap_1, y = umap_2),
            arrow = arrow(type = 'closed', length = unit(.1, 'inches')))

sobj |>
  AddMetaData(osteo_pt) |>
  FeaturePlot('Lineage2', cols = c('lightgrey','red')) +
  ggtitle('Pseudotime: Lineage2 in Prrx1+ cells') +
  geom_path(data = sds@curves$Lineage2$s, aes(x = umap_1, y = umap_2),
            arrow = arrow(type = 'closed', length = unit(.1, 'inches')))

# LPS receptor -------
sobj_mo |>
  DotPlot(c('TLR4','MYD88','TIRAP','TRAF6','TRPM2'), cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Monocyte clusters', title = 'Expression of LPS receptor & TRPM2')

sobj_mo |>
  DotPlot(c('SPI1','CEBPA','CTCF','TRPM2'), cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Monocyte clusters', title = 'Expression of TRPM2 & its TF')

# NFKB ---------
sobj_mo |>
  filter(group == '4 hpi') |>
  DotPlot(c('NFKB1','NFKB2','TRPM2'))

id_type <- sobj_mo |>
  distinct(seurat_clusters, trpm2_type) |>
  mutate(id = seurat_clusters, subtype = trpm2_type, .keep = 'none')

pbmc_nfkb <- last_plot() |>
  pluck('data')

pbmc_nfkb |>
  as_tibble() |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  left_join(id_type) |>
  ggplot(aes(TRPM2, NFKB1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'BM monocyte clusters 4hpi LPS',
       y = 'Average expression of NFKB1', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('lps.BM.mono.trpm2.nfkb1.cor', width = 70)

pbmc_nfkb |>
  as_tibble() |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  left_join(id_type) |>
  ggplot(aes(TRPM2, NFKB2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'BM monocyte clusters 4hpi LPS',
       y = 'Average expression of NFKB2', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('lps.BM.mono.trpm2.nfkb2.cor', width = 70)

